Numerical studies of confined states in rotated bilayers of graphene 
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Rotated graphene multilayers form a new class of graphene related systems with electronic prop- 
erties that drastically depend on the rotation angles. It has been shown that bilayers behave like 
two isolated graphene planes for large rotation angles. For smaller angles, states in the Dirac cones 
belonging to the two layers interact resulting in the appearance of two van Hove singularities. States 
become localised as the rotation angle decreases and the two van Hove singularities merge into one 
peak at the Dirac energy. Here we go further and consider bilayers with very small rotation angles. 
In this case, well defined regions of AA stacking exist in the bilayer supercell and we show that 
states are confined in these regions for energies in the [—7*, +7*] range with 74 the interplane mean 
interaction. As a consequence, the local densities of states show discrete peaks for energies different 
from the Dirac energy. 

PACS numbers: PACS number: 73.20.-r, 73.20.At, 73. 21. Ac, 61.48.Gh 



I. INTRODUCTION 

Graphene exceptional electronic properties are pre- 
dicted for an isolated layer of C atoms arranged on a 
pristine honeycomb latticejii^ Stacking several layers on 
top of each other, may affect the original electronic struc- 
ture. Indeed AB or Bcrnal stacking (such as what is 
found in graphite) destroys both the linear dispersion 
and the chirality properties even in a bilayer <^"— Many 
experimental realisations of graphene lead to the forma- 
tion of multilayers (on SiG^"— but also on metal surfaces 
such as Nii^ and in exfoliated flakesi^) where interac- 
tions between successive layers play a crucial role. While 
on the Si face of SiC, multilayers have an AB stacking 
and do not show graphene properties^'"— on the C-face, 
multilayers have been shown to present graphene like 
properties even when they involve a large number of C 
planes. ARPES^^i^ STM^ transport^ and optical 
transitionsi^ indeed show properties characteristic of a 
linear graphene like dispersion. These multilayers are ro- 
tated with respect to each other and the rotations show 
up as moire patterns on STM imagesJ^^i^ This ap- 
parent controversy -thick multilayers exhibiting graphene 
properties- was partially solved recently when different 
theoretical approachesi^^"— showed that rotated mul- 
tilayers are decoupled, at least for large rotation angles. 
Going further, theory predicts the existence of three do- 
mains that in fact correspond to two regimes: for large 
rotation angles, > 15°, the layers are decoupled and 
behave as a collection of isolated graphene layers. For 
intermediate angles, 2° < 6 < 15°, the dispersion re- 
mains linear but the velocity is renormalised. For these 
large and intermediate angles, Dirac cones persist and 
the interaction between rotated layers is a perturbation. 
What happens at the smallest values of 6 is more puz- 
zling. As already shown by different theoretical groups, 
for the lowest 6, flat bands appear and result in electronic 
localization. 



From the experimental point of view, while some ex- 
periments with graphene on SiG show none of the effects 
predicted by theory (no renormalisation nor van Hove 
singularities) r^ others^ do find them. Furthermore, 
Landau Level (LL) Scanning Tunneling Spectroscopy 
(STS) gave results in close agreement with theory and the 
three domains were observed for GVD graphene grown on 
Ni^ 

Here we go further and show that in the very small 
angle regime, confinement occurs and states become con- 
fined in the AA (all atoms on top of each other) region. 
Furthermore, confinement effects lead to the formation 
of sharp peaks in the local density of states in an energy 
window that is fixed by the interplane mean interaction. 
Thus, this apparently simple system -a graphene bilayer- 
presents a complex behavior ranging from ballistic to con- 
fined electrons when the rotation angle is changed. 

Tackling very small rotation angles - smaller than 1° - 
means handling very large cells that involve a huge num- 
ber of atoms - more than 60 000 -. This cannot be done 
from ab initio calculations and we developped a tight 
binding scheme^i using only p^ orbitals to do so. Ro- 
tated multilayers are described in part II, the method 
used is outlined in the third part and described in more 
details in the Appendix A. Part IV gives the variation of 
the van Hove singularity position with the rotation angle 
and parts V and VI describes the electronic structure of 
graphene rotated bilayers in the case of very small rota- 
tions angles. 



II. ROTATIONALY STACKED 
COMMENSURATE BILAYERS 

In the following we consider two graphene layers ro- 
tated in the plane by an angle 9. We start from an AA 
bilayer and choose the rotation origin O at an atomic site. 
Since we want to perform ab initio and Tight Binding cal- 



culations, we need a periodic system. A commensurate 
structure can be defined if tlie rotation clianges a lattice 

vector OB {m, n) to OB' {n, m) witli n, m tfie coordi- 
nates witli respect to tlie basis vectors ai (a/3/2,— 1/2) 
and 0,2 (■\/3/2, 1/2). Tlic rotation angle is then defined 
as follows: 



Cone 1 Cone 2 



cos6' 



Anm 



2{n? + nm + im?) 
and the commensurate cell vectors correspond to 



(1) 



t = OB' ~ ndi + 77102 , t' = —mdi + (n + m)a2. (2) 

The commensurate unit cell contains N — A{n'^ + nm + 
m?) atoms. As we have already shown,— the rotation 
angle is a good parameter to describe the system but the 
number of atoms is not since cells of equivalent size can 
be found for different angles. Indeed, large cells can be 
obtained for ~ - large n and m and small |?7i — n| - 
but also for large angles ~ 30° - then jm — n| is large- 
ly A good way to name the rotated layers is to give the 
(n, m) couple defined above, which is what we use in the 
present paper. 6 = 60° is the perfect AB stacking. 6 
close to 60° is obtained for n = 1 (tti = 1) and large m 
(n), it is equivalent in our calculations to ~ 0°. 



III. METHODS 

Our approach combines ab initio and tight binding cal- 
culations. First principle approach can tackle cells with 
up to ~ 1000 atoms. The tight-binding (TB) scheme al- 
lows electronic structure calculations by diagonalisation 
in reciprocal space for crystal structures containing up to 
59 644 atoms in a unit cell - corresponding to 6* = 0.469° 
and the (70, 71) bilayer -. By the using recursion tech- 
nique in real space, very large structures up to several 
millions - 3 006004 for the (500,501) with 6 = 0.067° - of 
atoms are studied. The tight binding scheme is described 
in Ref. |3l| and in the Appendix; it is restricted to Pz or- 
bitals. It was developped from ab initio calculations and 
it gives a good description of the rotated bilayers, for 
small angles as well as large ones but also for trilayers, 
graphene and graphite. The TB parameters were fitted 
to reproduce the velocity of monolayer graphene and the 
intcrlayer mean interaction jt = 0.34 eV (see Appendix). 
Flat layers were used. We performed an ab initio calcu- 
lation for the (6,7) bilayer allowing all the atoms to relax 
and found only minor changes: the layer corrugation is 
small and the effect on the band structure is very small»^ 
Velocities are calculated from the slope of the dispersions 
at the Dirac point along FK. 



IV. VAN HOVE SINGULARITIES 

The rotated bilayer has already been the subject 
of many theoretical developments either analytical or 
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FIG. 1: (color on line) Van Hove singularities of a (4,5) bi- 
layer. (a) schematic drawing, (b) (4,5) band structure and 
total DOS. 



numericalJ^^^— The main results are only briefly re- 
called in the following since we want to focuss on the 
very small rotation case and the confinement effect that 
results. These theoretical predictions are in very good 
agreement with Scanning Tunneling Spectroscopy exper- 
iments on CVD graphene grown on nickelJ^ Properties 
are symmetric with respect to 6* = 30°: In our calcula- 
tions, the behavior is the same for angles decreasing to 
0° or increasing to 60°. 

For large angles, 9 > 15°, the rotated bilayer behaves 
like two isolated graphene layers. For intermediate an- 
gles, 2° < 6* < 15°; the dispersion remains linear but 
the velocity is strongly renormalised. In this regime, one 
important feature is the appearance of van Hove singu- 
larities on both sides of the Dirac energy. In the follow- 
ing, we consider that layer 2 is rotated by an angle 9 
with respect to layer 1. Layer 2 Brillouin zone is then 
rotated by the same 9 with respect to layer 1 Brillouin 
zone. The closest Dirac cones belonging to layer 1 and 2 
intersect at a point P which is at mid-distance between 
the Ki and K2 points of the 1x1 Brillouin zones (see 
Fig.[T]a,). P is backfolded to point M of the supercell bril- 
louin zone. Interaction of the states with similar energies 
opens two gaps at the saddle point M (Fig.[T|D). Flat 
bands have a drastic effect on the electronic structure of 
a two-dimensionnal systems and appear as van Hove sin- 
gularities on each side of the Dirac energy (Fig.[T)D). In 
fact, each gap creates two van Hove singularities: Sin- 
gularities closest to the Dirac point appear as peaks _B-|_ 
and _B_ in Fig.[T|3, they are followed by a dip and a sharp 
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FIG. 2: (color on line) Energies of the van Hove singularities 
(the two peaks closest to Ed) in {n,m) bilayers versus 9). 
Lines are guide for the eyes. 
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FIG. 3: Gap at P (Point M of the supercell) as a function of 
6. Line is a guide for the eyes. 



increase of the DOS that corresponds to the flat bands 
on the other side of the gaps. The gap width and then 
the energy difference between the steep increase in the 
DOS and the peak arc equal to 2te with tg the mean 
interaction between states related to the two layers at 
the intersection between the Dirac cones, as defined by 
Ref.H. 

In Fig. [21 the energies, -E+ and E- of the closest peaks 
to the Dirac energy are plotted versus 9. For intermedi- 
ated angle values, 2 < < 15°, £'+ and E- are linear 
with respect to 0. 

As the rotation angle decreases, the intersections of 
the two Dirac cones belonging to the two graphene layers 
come closer to the Dirac energy (Fig.[T]). The intersection 
point P is at mid-distance between Ki and A'2. K1K2 
distance decreases linearly with 0. Since in this regime 
the dispersion is linear and the gap at P is roughly con- 
stant (Fig. [3]) , the position of the van Hove singularities 
then varies linearly with 6. 

For small angle values, 6^2°, £'+ and E- merge and 
form a sharp peak at the Dirac energy £'d . 

As shown in our previous article,— this sharp peak 
corresponds to states strongly localised in AA stacking 
regions of the moire. Recently Lopes Dos Santos et ali^ 
extended their treatment based on the continuum model 



to smaller angles and confirmed the localisation of the 
states in the AA regions of the moire pattern. 



V. CONFINEMENT IN AA REGIONS 

It is interesting to look at the band structures of bi- 
layers with small angles -small enough to have van Hove 
singularities at the limit of the linear variation (Fig.|3|)-. 
(19,20) bilayer is at the edge of the linear regime: the 
band structure still shows a linear dispersion and van 
Hove singularities. For the (31,32) bilayer, bands are 
flat and velocity is equal to zero^ but in the case of the 
(45,46) bilayer, the peak at i?D is split in two and the 
velocity at the Dirac point increases again and is equal 
to ~ 20% of the graphene velocity. This is coherent with 
the increase of velocity between magic angles found by 
Bistritzer and MacDonaldS^ (section IVl)) . 

It is important to notice that near the Dirac point one 
observes the existence of two bands that disperse linearly. 
These bands have a twofold degeneracy (Dirac cones from 
the two layers). The 4 bands contributes to the central 
peak of the total density of states. The weight of the 
eigenstates close to the Dirac energy is strong in the AA 
zone as shown by diagonalizationi '^^i^^ 

The analysis of the density of states, shows that a spe- 
cific regime occurs for very small angles, typically 6 < 1°. 
The total DOS (Fig.[S^) contains sharp peaks at energies 
different from the Dirac energy. Figs.[SjD and [St show the 
local DOS at the center of an AA stacked zone for de- 
creasing 0, and Fig.[SJl the local DOS at the center of an 
AB stacked zone. 

For so small rotation angles, the moire pattern presents 
well defined regions of AA and AB stacking. In the cen- 
ter of these regions, the local DOS follows the AA or AB 
DOS (dotted lines in Figs.[S]3, c and d), except for en- 
ergies close to the Dirac energy -. Indeed according to 
Fig. [5] peaks appear in an energy window defined by the 
mean interplane interaction 74. DOS in intermediate re- 
gions (not shown) are very similar to DOS in AB regions 
and do not show large peak. 

These results are consistent with confinement effect 
that can explain the discrete character of the spectrum 
in the [—7*, -l-7t] range. Indeed, in large moire patterns, 
AA regions seem isolated by AB and intermediate re- 
gions. From Fig.[5ji, we can say that confinement in AB 
regions is much less efficient than in A A regions. In a 
nanostructure, confinement results in peak in the DOS 
whose position depends on the size of the nanostructure. 
Here the size of the AA region is inversely proportional 
to the rotation angle 0. In Fig. [Hi we plot the energy of 
the first peak away from the Dirac energy as a function 
oi 0. As expected this energy decreases when the size 
of the A A region increases, which is when the angle 
decreases. 

In the energy range [— 7t,+7t], the symmetry of the 
wavefunction is different in an AA or AB bilayer (Fig. [3 
see also the DOS in the appendix Fig.llOb). Because of 
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FIG. 4: (color on line) Band structures and DOS for bilayers 
with rotation angles that lead to van Hove singularities close 
to the limit of the linear variation, (a) (19,20), 6 — 1.67°, 
(b) (31,32), e = 1.05°, (c) (45,46), 9 = 0.73°. Left: Disper- 
sion of monolayer graphene is shown in black dashed line for 
comparison. Right: (red line) local DOS at the center of AA 
zone, (black dashed line) total DOS. 



the different symmetry of the wavefunctions, matching is 
difficuh which could explain that the AA stacked region 
is isolated. We note also that in this energy range the 
coupling between electrons and holes which have oppo- 
site velocities could favor localization. On the contrary 
out of the energy range [— 7t,+7t], the mixing between 
the hole states of one plane and the electron states of 
the other plane is weak. Indeed the coupling occurs be- 
tween electron states of one plane and hole states of the 
other plane that have the same wave vector but have a 
difference in energy that is much larger than 74. There- 
fore outside of the energy range [—7*, +Jt] the mixing of 
states of the two planes occurs separately between the 
electron states of the two planes and between the hole 
states of the two planes. We suggest that this absence of 
mixing between electrons and holes which have opposite 
velocities suppresses the localization. 
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FIG. 5: (color on line) TB DOS in bilayers: (a) Total DOS. 
Local DOS: (b) (c) on atom at the center of an AA zone, (d) 
on atom at the center of an AB zone (At. A: with an atom 
directly on top or below and at. B: without an atom on top or 
below). {n,m) bilayers: values (100,101) 9 = 0.33°, (200,201) 
9 = 0.16° and (300,301) 9 = 0.11°. (400,401) 9 = 0.08°. 



VI. MAGIC ANGLES 

As discussed previously, at very small angles, the ve- 
locity of rotated bilayers is no longer monotonic but 
shows minima and maxima (Fig. [8]). Minima correspond 
to a flat band with zero velocity at the Dirac energy. 
Bistritzer and MacDonald in Rcf. [23 found that the ve- 
locity is equal to zero for special magic angles 0„. The 
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FIG. 6: Variation of the energy of the first peak (-Bi) away 
from the Dirac point of the local DOS at the center of AA 
region as a function of 6. 



values in degrees are 9i — 1.05, 62 = 0.5, ^3 = 0.35, 
04 = 0.24, 6*5 = 0.2. We note that this series is simply 
given by 6'„ = 1.05/n. In the present calculations we 
recover by a full tight-binding calculation the first two 
angles n = 1 and n = 2 (see Fig. [8]) which are given by 
9n — 1.13/?7 in good accordance with the more simplified 
model^ of Ref.M 

Below we propose a simple heuristic argument which 
relates these magic angles to quantization conditions. At 
these magic angles the bands become flat which means 
that the states of the different AA zones are decoupled 
and thus confined in the AA zones. In the following, we 
use a hard wall model to describe the confinement effect. 
The confinement condition is typically: 
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fed + 2$ = 2n7r, 



(3) 



where k is the wavevector counted from the K point of 
the Dirac wave in the AA zone, | is the typical size of the 
AA zone and $ is the phase accumulated at the reflexion 
at the boundary of the AA zone. We assume that this 
phase can be integrated in an efficient size of the AA zone 
so that the quantization condition finally gives: 



xkD ~ 2mr, 



(4) 



FIG. 7: (color on line) Band structure of an AA and an AB 
stacked bilayer in the vicinity of the Dirac Energy. 
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FIG. 8: (color on line) Velocity of a bilayer divided by the 
velocity of a monolayer, at the Dirac point, as a function of 6. 
The insert shows a zoom on very small angles and on angles 
close to 60°. Black/green point: TB calculation, red cross: 
ab initio calculation, blue line: law proposed by Lopes dos 
Santos et al.,— orange line: model proposed by Bistritzer and 
MacDonalcl^S, for very small angles. The latter line is rescaled 
to obtain a first minimum of the velocity at 1.13° as in TB 
calculations (see text). Ab initio calculations and most of TB 
calculations for 9 > 1.05° are from Ref. [3ll. 



where D is the period of the moire and a; is a numerical 
coefficient of the order of 1 such that xD/2 is the effec- 
tive size of the AA zone. This expression can be cast in 
terms of the parameters of the present study. Indeed the 
period D of the moire expressed in nanometers and the 
angle 9 of the rotation expressed in degrees are related 
by £'(nm) = „,"', a for small angles. In addition, we have 
the relation hvpk = 74 between the parameter 74 and the 
wavevector k at the Dirac energy in the AA zone. Then 
the quantization relation leads to: 



3.4x74 (eV) 



^(deg). 



(5) 



For the present value of the parameter 7^ = 0.34 eV and 
with a reasonable value of the parameter x ^ 0.95 this re- 
lation leads exactly to the magic angles 9{deg) = 1.13/n 
found here. This value is close to the series 1.05/n found 
in Ref. {M- 

Why the simple hardwall argument discussed here can 
describe the magic angles series is not obvious. Indeed 
Bistritzer and MacDonald^S used an interplane interac- 
tion hamiltonian with a sinusoidal variation in space. 
For free electrons a sinusoidal interaction potential would 
probably not lead to the localization with hard wall con- 
ditions. Yet as discussed above, the localization in the 
energy range [—74, +^t\ that is found in this work suggest 
that the spinor character of the electronic wavefunction 
in graphene leads to more complex localization phenom- 
ena. 



VII. CONCLUSION 

Thanks to the tight binding scheme that we devel- 
opped, we studied the electronic structure of rotated bi- 
layers in the case of very small rotation angles, which 
is beyond what was recently proposed. We showed that 
states are confined in regions of well defined AA stacking 
and as a consequence that the LDOS shows discrete peaks 
not only at the Dirac energy but in the [—74, +7t] range 
with 74 the mean interlayer hopping integral. Very large 
nioire patterns (corresponding to angles close to 1°) have 
already been observed experimentally. This system could 
then offer the possibility to study the physics of a regular 
array of weakly interacting dots. To conclude we note 
that the rotated bilayers of graphene present analogies 
with quasicrystals and approximant crystals that have 
also large unit cells and tend to confine electrons (Ref. 
I37I and references within) . In the context of quasicrystals 
this localization by large unit cells leads to unique trans- 
port properties.—!^ This suggests that rotated graphene 
bilayers with large unit cells could have original transport 
properties. 
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Appendix A: Tight Binding scheme for bilayers 

In this appendix we present in details the tight binding 
(TB) scheme that we have used in our previous work^ 
and in the present results. It reproduces the ab initio 
calculations of the electronic states for energies within 
±2 eV oi Et) whatever the rotation angle is (Fig. 13 [10] 
and Ref. \3u \. Here Ed is set to zero. 

Only pz orbitals are taken into account since we are 
interested in what happens at the Fermi level. Since the 
planes are rotated, neighbours are not on top of each 
other (as it is the case in the Bernal AB stacking). In- 
terlayer interactions are then not restricted to ppa terms 
but some ppn terms have also to be introduced. The 
Hamiltonian has the form : 



H = 5I^^N)(*i + E ^'jN)01' 
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FIG. 9: (color on line) Band structure E{k) calculated from 
ab initio VASP and TB: (a) graphene, (b) AB bilayer, (c) AA 
bilayer. 74 — 0.34 eV. 




where \i) is the Pz orbital located at r^, and {i,j) is the 
sum on index i and j with i ^ j. The coupling matrix 



FIG. 10: (color on line) Density of states: (a) ab initio (Linear 
Muffin Tin Orbitals^ and TB total DOS in graphene, (b) TB 
total DOS in AB (Bernal) and AA bilayers. [Insert: total 
DOS arround the j5d = 0. 7* = 0.34 eV] 
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FIG. 11: Coupling term t between two pz orbitals versus 
the distance r: (full line) in two layers (inter-layer coupling); 
(dashed line) in the same layer (intra-layer coupling). Here, 
the smooth cut-off function Fc is not taken into account (ie 
large Tc value in (|A7p ). 



where a is the nearest neighbor distance within a layer, 
a = 1.418 A, and ai is the interlayer distance, ai = 
3.349 A. First neighbors interaction in a plane is taken 
equal to the commonly used value, 70 = 2.7 eVi^ Second 
neighbors interaction 7q in a plane is seii^ to 0.1 x 70 which 
fixes the value of the ratio q-n/a in equation ()A4|) . The 
inter-layer coupling between two pz orbital in tt configu- 
ration is 71 . 71 is fixed to obtain a good fit with ab-initio 
calculation around Dirac energy in AA stacking and AB 
bernal stacking and then to get 74 = 0.34eV (figures[71[n] 
and fTD)) which results in 71 = 0.48 eV. We chose the same 
coefficient of the exponential decay for V^j,^ and Vpp^, 
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= 2.218. 
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element, iy , between two Pz orbitals located at fi and r. 
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40 



Uj = n^Vppairij) + {l-n^)Vpp.„{rij), (A2) 



where n is the direction cosine of r^ = 

axis and Tjj the distance r^ between the orbitals. 



-^ and r 



^3 



■ fi along (Oz) 



(A3) 



Zij is the coordinate of fij along (Oz). It is either equal 
to zero or to a constant because the two graphene layers 
have been kept flat in our model4i We use the same fol- 
lowing dependance on distance of the Slater and Koster 
parameters: 



yppTT vij ) ' 
Vppairij) = 71 e' 



,9. 1- 



70 e 
9^(1 



->F,{n,), (A4) 

FciUj)- (A5) 



with flo = 2.456 A the distance between second neighbors 
in a plane. Intra-layer and inter-layer coupling terms 
t{rij) calculated with these parameters are shown on fig- 
ure [TTJ In (jA4p and (|A5|) . a smooth cut-off function is 
introduced)^ 



Fcir) = [l + e^^j , (A7) 

with Tc the cut-off distance and Ic = 0.5 a.u. = 0.265 A.— 
For r <C Vc, Fc{r) ~ 1; and for r ^ Tc, Fc{r) ~ 0. 
All results presented in this article are calculated with 
Tc = 2.5ao = 6. 14 A. We have checked that the results 
are independent of r^ for r,, > 4.9 A. 

All Pz orbitals have the same on-site energy e^ (equa- 
tion (|A1[) ). ei is set to —0.78 eV so that the energy _Ed of 
the Dirac point is equal to zero in mono-layer graphene. 
ei is not zero because the intra-layer coupling between 
atoms beyond first neighbors breaks the electron/hole 
symmetry and then shifts Ed. 
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